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ABSTRACT 

Gravitationally unstable accretion disks emerge in a variety of astrophysical contexts — giant planet 
formation, FU Orioni outbursts, feeding of AGNs, and the origin of Pop III stars. When a gravita¬ 
tionally unstable disk is unable to cool rapidly it settles into a quasi-stationary, fluctuating gravito- 
turbulent state, in which its Toomre Q remains close to a constant value Q o ~ 1. Here we develop 
an analytical formalism describing the evolution of such a disk, which is based on the assumptions of 
Q = Q o and local thermal equilibrium. Our approach works in the presence of additional sources of 
angular momentum transport (e.g. MRI), as well as external irradiation. Thermal balance dictates a 
unique value of the gravitoturbulent stress a g t driving disk evolution, which is a function of the local 
surface density and angular frequency. We compare this approach with other commonly used gravito¬ 
turbulent viscosity prescriptions, which specify the explicit dependence of stress a g t on Toomre Q in 
an ad hoc fashion, and identify the ones that provide consistent results. We nevertheless argue that 
our Q = Q 0 approach is more flexible, robust, and straightforward, and should be given preference 
in applications. We illustrate this with a couple of analytical calculations — locations of the snow 
line and of the outer edge of the dead zone in a gravitoturbulent protoplanetary disk — which clearly 
show the simplicity and versatility of the Q = Q o approach. 

Subject headings: accretion, accretion disks — instabilities — (stars:) planetary systems: protoplan¬ 
etary disks — (galaxies:) quasars: general 


1. INTRODUCTION. 

Astrophysical accretion disks can be prone to gravi¬ 
tational instability (hereafter GI) when their tempera¬ 
ture is low and surface density is high (Safronov 1960; 
Toomre 1964). In particular, GI is conceivable in the 
outer parts of protoplanetary disks (Cameron 1978; Boss 
1998) at the early stages, when the disks are still mas¬ 
sive because of ongoing infall. Accretion outbursts of FU 
Orioni stars can be driven by the gravito-magnetic cycle 
in the dense, gravitationally unstable parts of the proto¬ 
planetary disk (Audard et al. 2014). Outer parts of AGN 
disks accreting at high M are also expected to become 
gravitationally unstable far from the black hole (Paczyn- 
ski 1978a, 1978b; Kozlowski et al. 1979; Kolykhalov & 
Sunyaev 1980; Goodman 2003). Turk et al. (2009), Clark 
et al. (2011), Greif et al. (2012) find gravitationally un¬ 
stable disks around young Population III stars in their 
simulations of star formation at redshifts of z ~ 20 — 30. 

Possibility of the GI in a gaseous disk is characterized 
by the so-called Toomre Q parameter defined as (Toomre 
1964) 


where E and c s = (fcsT/V) 1 / 2 are the surface density 
and sound speed of the disk, and fl = ( GM/r 3 ) 1 / 2 is the 
angular frequency (assuming Keplerian rotation around 
a central object with mass M). In the linear regime GI 
sets in as Q —> Q o from above, where Q o ss 1 — 1.5 is 
the threshold value suggested by numerical experiments 
(Boss 2002; Cossins et al. 2009, 2010). 

Operation of the GI is accompanied by enhanced angu¬ 
lar momentum transport in the disk driven by the non- 
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axisymmetric gravitational torques. As a result, GI re¬ 
sults in mass redistribution in the disk, driving its evo¬ 
lution. Energy dissipation caused by the application of 
gravitational stresses heats up the disk and tends to op¬ 
pose the GI. 

Because of this feedback, the nonlinear outcome of the 
GI sensitively depends on another dimensionless param¬ 
eter — the product of the local cooling time in the disk 
Cool and fl. Gammie (2001) demonstrated that when 
cooling is fast and flt coo \ < /3 cr it ~ 1, the disk disinte¬ 
grates into a number of bound self-gravitating structures. 
Such fragmentation has been invoked by Boss (1998) to 
explain the origin of giant planets (see Rafikov (2005, 
2007) for constraints on this scenario). Goodman & Tan 
(2004) and Levin (2007) suggested that fragmentation 
of gravitationally unstable quasar disks should result in 
formation of massive stars migrating through the disk. 
Clark et al. (2011), Greif et al. (2012), Latif & Schleicher 
(2014) propose that fragmentation of protostellar disks 
around Pop III stars can produce low mass extremely 
metal poor stars that can survive until present days. 

In the opposite limit of long cooling time flt coo \ ^ /3 C rit 
it was shown by Gammie (2001) that the disk settles 
into a quasi-stationary state of gravitoturbulence. In this 
regime surface density fluctuates in time, but the disk 
state averaged over the period longer than the dynami¬ 
cal timescale U” 1 remains the same. The time-averaged 
value of the Toomre Q parameter in a gravitoturbulent 
disk is around Q o, i.e. the disk maintains itself in a 
marginally stable state. This general picture has been 
confirmed with global simulations by different groups 
(Rice et al. 2003, 2005; Durisen et al. 2007; Cossins et 
al. 2009,2010). 

The critical value of the cooling time /3 cr it^ _1 corre¬ 
sponding to the transition between the gravitoturbulent 
and fragmenting regimes, depends on a variety of fac- 
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tors. One of the key determinants is the equation of state 
(EOS) of the gas, with softer EOS promoting fragmenta¬ 
tion and resulting in higher /3 cr it (Rice et al. 2005; Jiang 
& Goodman 2011). Opacity transitions, e.g. due to dust 
grain evaporation, also affect critical t coo i (Johnson & 
Gamrnie 2003), as well as other forms of the tempera¬ 
ture dependence of opacity (Cossins et al. 2010). Exter¬ 
nal irradiation (Rice et al. 2011) and the details of the 
disk structure (Meru & Bate 2011a) may also affects the 
value of /3 C rit- 

Some concerns have been raised regarding the conver¬ 
gence of gravitoturbulent disk simulations and the ex¬ 
istence of the well-defined /3 cr it, as its value has been 
claimed (Meru & Bate 2011b; Paardekooper 2012) to 
vary with the grid resolution. However, later this non¬ 
convergence has been traced mainly to the numerical ef¬ 
fects (Paardekooper 2011; Meru & Bate 2012; Rice et al. 
2014). 

It has also been debated whether fragmentation ensues 
because rapid cooling facilitates collapse of unstable frag¬ 
ments, or because the disk can withstand only certain 
amount of stress and fragments when the dimensionless 
stress parameter a (Shakura & Sunyaev 1973) exceeds a 
threshold value a C rit • Since in thermal equilibrium in the 
absence of external energy inputs (Gamrnie 2001) 

a = 9 7 ( 7 - 1) (^ c ° o1 ) ’ 

where 7 is the adiabatic index of gas, the values of a cr it 
and /3 cr it are directly connected. Comparing simulations 
with different 7 (different EOS) Rice et al. (2005) noted 
that a cr it is essentially independent of 7 and is about 
0.06. This led them to suggest that the primary reason 
for fragmentation is the maximum stress that can be sus¬ 
tained by the disk. However, simulations with external 
irradiation (Rice et al. 2011) suggest that a cr i t does de¬ 
pend on the level of irradiation and is thus non-universal. 

The main goal of this work is to explore disk charac¬ 
teristics in the gravitoturbulent state, which can persist 
for a long time over a significant range of radii. For that 
reason it is important to understand how the disk evolves 
in this regime under the action of the non-axisymmetric 
gravitational torques. Currently direct 2D and 3D sim¬ 
ulations are too numerically expensive to permit such 
exploration, and one often has to resort to azimuthally- 
averaged, one-dimensional (in radius) disk models. To 
evolve them properly one must provide a description of 
the angular momentum transport by the gravitoturbu- 
lence. Formulating such a description is the focus of the 
present work. 

Balbus & Papaloizou (1999) argued that due to the 
long-range nature of the gravitational interaction the an¬ 
gular momentum transport due to the non-axisymmetric 
gravitational perturbations is inherently non-local. How¬ 
ever, Gamrnie (2001) showed that in cold, geometrically 
thin disks angular momentum transport by the gravita¬ 
tional torques can still be described as a local process. 
Lodato & Rice (2004, 2005) and Cossins et al. (2009) 
confirmed this conclusion numerically, certainly for the 
low-mass disks. 

In this work we adopt the latter point of view and will 
explicitly assume the gravitoturbulent transport to be 
local and characterized by the effective viscosity param¬ 
eter a g t- Different explicit and implicit prescriptions for 


a g t have been proposed, which can be classified into two 
general categories. 

The first class of a g t prescriptions relies on the fact 
that in local dynamical and thermal equilibrium the an¬ 
gular momentum transport in the disk is intimately re¬ 
lated to its thermal state. This allows one to directly 
express a g t via the disk temperature and density. Using 
constant M disk without external energy inputs as an ex¬ 
ample, Rafikov (2009) has demonstrated that this prop¬ 
erty, when coupled with the defining characteristic of the 
gravitoturbulence — the condition Q « Qq, allows one 
to directly relate a g t to E providing a complete descrip¬ 
tion of the disk evolution. This approach has also been 
implicitly featured in several numerical studies (Terquem 
2008; Clarke 2009; Zhu et al. 2009a). 

A different way of describing the gravitoturbulent disk 
evolution does not explicitly assume Q ss Qq. Instead, it 
specifies the explicit dependence of a g t on Q (Kratter et 
al. 2008; Zhu et al. 2009b, 2010a, 2010b; Martin & Lubow 
2011, 2013, 2014). This approach dates back to the work 
of Lin & Pringle (1987), who suggested that a gt ~ Q~ 2 
based on a phenomenological description of transport in 
gravitationally unstable disks. However, in most cases 
such prescriptions do not naturally follow from physical 
arguments. Instead, they are designed to replicate cer¬ 
tain qualitative features of a g t behavior. 

The goal of this work is to formulate, analyze and com¬ 
pare different approaches to characterizing a gt . First, 
in we discuss equations used to describe evolution 
of gravitationally unstable disks. Then, in f[3] we high¬ 
light the differences between the two aforementioned ap¬ 
proaches to closing the system of evolution equations. 
After considering the steady, constant M models in JU 
we contrast the performance of different gravitoturbu¬ 
lent viscosity prescriptions in fj5) We discuss our results 
in (j 6 j where we show, in particular, how the gravitotur¬ 
bulent prescription based on Q ss Qq condition can be 
naturally used to obtain simple analytical expressions for 
the locations of the ice lines and dead zone edges in the 
gravitoturbulent protoplanetary disks. 


2. BASIC EQUATIONS. 


Provided that gravitoturbulent stress can be described 
as local a-viscosity, it is convenient (Rafikov 2013) to 
characterize the disk structure via the angular momen¬ 
tum flux Fj defined (in a Keplerian disk) as 


Fj = 37 tuEZ = 37 raCgE r 2 , (3) 

where l = fir 2 is the specific angular momentum and v 
is the kinematic viscosity expressed through the dimen¬ 
sionless parameter a as (Shakura & Sunyaev 1973) 



As always, we will assume the gas sound speed c s to 
be determined by the midplane temperature of the disk. 
Angular momentum transport is effected both by grav¬ 
itoturbulence and by any other background sources of 
effective viscosity, such as MRI. 

Mass accretion rate through the disk M (defined to be 
positive for inflow) is related to Fj via a simple relation 
(Rafikov 2013) 


M = 


dFj 
dl ' 


(5) 
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Coupled with the continuity equation, this results in the 
following evolution equation: 

<9£ _ 1 d 

dt 2nr dr 

Expressing Fj via equation © one reproduces the clas¬ 
sical equation for the viscous disk evolution (Lightman 
& Eardley 1974; Lin & Papaloizou 1996) 


(m\ 

1 dFj 

\dr J 

dr 


( 6 ) 


as 

~dt 


3d_ 
r dr 



(7) 


Solving either of these equations requires the knowledge 
of the thermodynamical properties of the disk since both 
Fj and v are functions of gas temperature T, see equa¬ 
tions © & ©• 

We now address the thermal state of the disk. We 
assume the disk to be in thermal equilibrium, so that en¬ 
ergy sources and sinks are in local balance. The former 
include viscous dissipation and external irradiation in¬ 
tercepted by the disk. The latter is effected by radiative 
cooling. For local angular momentum transport the rate 
of viscous energy dissipation per unit radial distance in 
the disk is dE/dr = —FjdVt/dr. Then the rate of energy 
production by viscous stresses per unit area of the disk 
is 


1 dE _ 3 Fjfl 
2irr dr 47 r r 2 


( 8 ) 


Gravitationally unstable disks are often studied in the 
limit of a self-luminous disk (Gannnie 2001; Rafikov 
2009), in which viscous heating of any nature dominates 
over the energy input due to external irradiation. This 
assumption should be reasonable for e.g. Class 0 T Tauri 
stars, which are characterized by intense mass accretion 
and relatively low luminosity of a very young central star, 
which is also heavily obscured. 

However, in general, the disk, especially its outer parts, 
may also be heated appreciably by the external radiation 
field with temperature T; rr (r) (Rice et al. 2011; Zhu et 
al. 2012); Ti rr can be a function of r if irradiation is 
due to the central object. To account for this possibility 
we describe the local thermal balance, which ultimately 
determines the midplane disk temperature T, via the fol¬ 
lowing approximate relation: 


Q = 


3 Fjfl 

4tT j-2 


2a 

/(t) 


(T 4 - 7j 4 r ) • 


(9) 


Details of the vertical transport of radiation in the disk 
are specified here via an explicit function f(r) of the 
optical depth 


t « k(T)£. (10) 

To zeroth order the disk optical depth is determined pre¬ 
dominantly by the midplane value of the temperature T 
(k can also be a function of gas density). 

The explicit form of /(t) depends on whether the disk 
is optically thick or thin and on the mechanism of the ver¬ 
tical energy transport (Rafikov 2007). Here we assume 
radiative energy transfer (Rafikov 2007 has derived /(t) 
for convective transport of energy) with dust-dominated 
opacity. In this case /(r) can be reasonably well approx¬ 
imated (up to constant factors of order unity) by 

/(t)~t + t _1 . (11) 


This expression does not discriminate between the Rosse- 
land and Planck mean opacities and is approximately 
valid (up to constant factors of order unity) both in the 
optically thick and thin regimes (Goodman 2003). It 
is easy to see that the condition © supplemented with 
equation CP properly describes the disk thermal balance 
in different asymptotic regimes of r. Indeed, in the opti¬ 
cally thick case (r 1) one finds aT 4 « erT( 4 r + (Q/2)t, 
while in the optically thin regime (r <C 1 ) one has 
aT 4 « erX| 4 r + (Q/2)t~ 1 . Both expressions are valid up 
to constant factors in the appropriate limits. 

For simplicity in this work we will consider the follow¬ 
ing behavior of k appropriate for dust opacity in certain 
temperature regimes (Bell & Lin 1994): 

k = k 0 T^, ( 12 ) 

with Ko; d being constants. For example, at low tem¬ 
peratures, below 170 K the opacity is dominated by ice 
grains and is characterized by (Bell & Lin 1994) 

(3 = 2, ko = 5x 10 -4 cm 2 g _1 K -2 . (13) 

Our results can be trivially extended to other, more com¬ 
plicated forms of n. 

Equations CT- (fl3l) provide the sought relation between 
the disk temperature and surface density, which is used 
to determine the viscosity behavior in equation Cl- 

Our final note on disk thermodynamics concerns the 
case when cooling is not due to dust emission. In particu¬ 
lar, Latif & Schleicher (2014a,b) considered the structure 
of gravitationally unstable disks around Pop III stars. 
Such disks are believed to be metal-free (i.e. no dust 
opacity), with cooling provided by molecular hydrogen. 
In this case, the right-hand side of the energy balance 
equation © should be modified to 2(c s /H)A, where c s /H 
is the disk scale height and A is the volumetric cool¬ 
ing rate. This expression is valid in the optically thin 
regime, with certain modifications required when line 
cooling switches to the optically thick regime (Latif & 
Schleicher 2014a; Ripamonti & Abel 2004). 

3. CLOSURE IN THE GRAVITOTURBULENT STATE. 

Definition © and thermal balance condition © sup¬ 
plemented with equations m and m allow one to 
uniquely relate T to £ if the behavior of a is known. 
This situation is typical for the disk regions dominated 
by the background viscosity (e.g. due to MRI), where one 
usually knows (or, more often, postulates) some behav¬ 
ior of a. Knowing the T(£) dependence one can express 
Fj and v as a function of £ only, leaving £ as the only 
dependent variable in the equation ©■ This provides a 
closure of the system of equations for the density, angu¬ 
lar momentum, and thermal balance, and allows one to 
self-consistently evolve £ in time and space. 

In gravitoturbulent state the behavior of a is not 
known a priori and this approach does not apply. To 
close the system of evolution equations one needs to make 
additional assumptions. Two conceptually different ap¬ 
proaches to this problem are described next. 

3.1. Closure via the Q = Qo condition. 

Rafikov (2009) pointed out that a necessary closure is 
naturally provided by the the basic property of a gravito¬ 
turbulent disk following directly from simulations — that 
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the disk hovers on the margin of gravitational stability 
with 


Q = Qo, (14) 

where Qo ~ 1- 

Indeed, with definition © this condition predicts an 
explicit relation between T and E in the form 


r Q (f,s) 


jj_ f nGQ 0 \ 2 2 

k B V H( r ) ) 


(15) 


Plugging this expression into the thermal balance con¬ 
dition © gives the following expression for the gravito- 
turbulent stress 




Stt r 2 a{T%-Tr) 
3 B(r) /(r) 


(16) 


The dependence on surface density E arises here because 
Tq, which enters this expression both explicitly and also 
through 


r(r,E) = E K (T Q (r,E)), (17) 


is a function of E. 

Combining equations © - (UGH one also comes up with 
the explicit expression for the effective a-parameter due 
to gravitoturbulence: 


«gt(cS) 


8 <t(7tGQo) 6 

9 /(r) 



The dependence on r comes only through fl (r). This 
result was first obtained by Rafikov (2009) in the non- 
irradiated case (Ti rr = 0) when studying constant M 
gravitoturbulent disks. However, it clearly also holds 
more generally, for evolving and irradiated gravitotur¬ 
bulent disks, because thermal balance gets established 
faster than the viscous evolution proceeds. Note that 
Rice et al. (2011) found a different (linear) scaling with 2 
Ti rr in their expression for the gravitoturbulent a be¬ 
cause they used a different cooling model, namely as¬ 
suming a constant cooling time. Expression (fl8l) for a^ t 
completes the closure and makes evolution equation © 
self-consistent. 

On the other hand, gravitoturbulent disk evolution can 
be described entirely without introducing a-parameter. 
Indeed, substituting Fj = Fj igt (r, E) given by the ex¬ 
pression (ITUD into equation © one finds 


dT, 

dt 


8 d 
3 r dr 


J_d_ 

Qr dr 


f 2 ° {Fq ~ ^irr) 1 \ 
ft f(r(T Q )) Jj’ 


(19) 


where, again, Tq = Tq{t, E), and T m {r) behavior is 
specified. This (in general nonlinear) equation is a 
closed-form, fully self-contained evolution equation for 
E as a function of t and r. It represents one of the main 
results of this work. 

In the parts of the disk where external irradiation 
plays the dominant role this equation adopts a simple 
time-independent solution Tq —>• Ti rr , so that E (r) = 
H(r)c s (Ti rr (r)) /(ttGQq) (Rafikov 2009). 


2 Note that the expression in parentheses in equation m can 
be written as (l — Qf rr /<9o), with Q; rr = Qe s (T; rr ) /(irGE) (Rice 
et al. 2011). 


3.1.1. Additional sources of viscosity. 

Angular momentum transport in the disk may be ef¬ 
fected not only by gravitoturbulence but also by addi¬ 
tional stresses. This would be the case, for example, if 
the disk is both sufficiently ionized for the MRI to oper¬ 
ate and massive enough for being gravitationally unsta¬ 
ble. 

To account for the possibility of additional, non- 
gravitoturbulent viscosity parametrized by a-parameter 
a m , one can simply write 

Fj =ipFj tg b(r, E) + Fj, m (r, E, T). (20) 

Here Fj gt (r , E) is still given by the expression (fTUD with 
the switch function ip introduced as described below. 
The non-gravitational viscous angular momentum flux 
is 

Fj,m{r, E, T) = 2>na m c? s (T)Yir 2 , (21) 

and the behavior of a m is specified. Thermal balance 
relation m now gives 

LT? ^ ^ m _ 8 tt r 2 a (T 4 - T^) ^ 

V’ J, gt(r, )+ j,m(r, , ) 3 n y( r ( E)T )) ( ) 

from which we determine T as a function of r and E. 
Then equation (BUI) yields Fj as a function of r and E, 
allowing the evolution of E to be followed using equation 

©■ 

A subtle point in this prescription is that equation 
(l 20 l) always includes the gravitoturbulent transport in 
the form m, with Tj igt independent of the actual disk 
temperature T and the value of Q. As a result, our pre¬ 
scription formally has non-zero gravitoturbulent stress 
Fj tgt even in the gravitationally stable part of the disk, 
which is not expected. 

For this reason we introduce a switch function il> in 
equation (l20l) . which takes care of this issue. There are 
several different ways in which this can be done. One 
method would be to adopt ip which quickly goes to zero 
as soon as Q > Qo- However, in the absence of a micro¬ 
scopic theory of gravitoturbulence such a factor would 
necessarily be introduced in an ad hoc fashion. 

Instead, we have chosen to use 

= 0(F J>gt ), (23) 

where 9(z) is the Heaviside step function (9(z) = 1 for 
z > 0; 9{z) = 0 for 2 < 0). With this approach we 
simply keep Tj, g t in the form (fTUl) in equation (l 20 l) . as 
long as it is positive. This is not a problem, since, as 
demonstrated in Rafikov (2009), as soon as Q exceeds 
the threshold value Qo, the transport is guaranteed to be 
dominated by the viscous stress, i.e. Fj tTn Fj tgt . In 
other words, a^ t <C a m in gravitationally stable disks, 
even if Ti rr = 0. And vice versa, Fj tTn -C Tj igt and 
a m when the disk is gravitationally unstable and 
Q —> Qo- This situation is further discussed in (©and is 
illustrated in Figure[5j where we display the actual run of 
the viscous and gravitoturbulent transport contributions 
for a particular disk model. 

Under certain circumstances one may find that Tq < 
Ti rr in the gravitationally stable part of the disk, so that 
Fj, g t < 0 according to equation (fTUl) . In this case our 
choice BUD of ip simply guarantees that the gravitotur¬ 
bulent part of the angular momentum flux vanishes and 
a = a m exactly. 
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Evolution equation © based on the prescription (TSUI) 
with %p given by equation (l23l) allows us to explore the 
transition between the disk regions dominated by gravi- 
toturbulence and background viscosity. In the appropri¬ 
ate limits its solution reduces to either the gravitoturbu- 
lent disk solution (for large r) following from equation 
(HU). or the conventional viscous disk solution (for small 
r, where Fj <C Fj^ m ). Alternatively, one can evolve 
the disk structure using equation © with v determined 

by 

a = ipo^ t (r, E) + a m . (24) 

It is easy to show that this approach results in almost 
the same disk properties, except for the region where 
Q!g t ~ a m and the disk is close to marginal gravitational 
stability. 

In summary, the prescription (12011 with Fj gt given 
by (flfil) should be applicable for any value of Q (even 
though only approximately in the parts of the disk just 
transitioning to the marginally gravitationally unstable 
regime). 


3.2. Closure via the explicit a g t(Q) dependence. 

An alternative closure scheme that has been broadly 
used in the literature postulates some dependence of the 
gravitoturbulent viscosity a g t on Q. In this case closure 
is analogous to a regular viscous disk anzatz: equations 
©, © with a = a m +a gt (Q), and (©-(HU are combined 
to yield a unique T(r, E) relation. It is then plugged into 
the relation (©]) for viscosity, resulting in the expression 
for v(r 1 E) and allowing equation © to be evolved in 
time. 

Note that the condition (1141) is not used explicitly in 
this approach and thus in general there is no guarantee 
that this prescription would result in a truly gravitotur¬ 
bulent disk structure, with the disk hovering at the edge 
of instability with Q Qq. 

In the absence of a microscopic model of gravitoturbu- 
lence that could motivate a possible dependence of a g t 
on Q , several ad hoc prescriptions for a g t(Q) have been 
suggested. Their main feature is the rapid increase of 
a g t as Q approaches some critical value ~ 1 from above. 

In particular, Zhu et al. (2010a,b) used 3 


«gt(Q) = e q4 


Martin & Lubow (2011,2013,2014) adopted 

i/ii \ 2 

a et(Q) = m & x { a. 


Q 


crit 


Q 


- 1 


,o 


(25) 


(26) 


with a m = 10~ 2 and Q C rit = 2. The Q~ 2 dependence 
is motivated by the work of Lin & Pringle (1987). Note 
that this prescription explicitly relates the value of grav¬ 
itoturbulent viscosity to the background viscosity a m . 
This recipe was also used in Owen & Jacquet (2014) to 
explore the chemical evolution of a protoplanetary disk 
undergoing accretion outbursts. 

Kratter et al. (2008) have used 


«g t (Q) = max < 0.14 


1.3 


max(<3,1) 


- 1 


,0 .(27) 


3 Zhu et al. (2009b) used a very similar prescription a g t(Q) = 
e-« 2 . 



Q 


Fig. 1.- Illustration of different explicit a g t(Q) prescriptions 
given by equations (l25l-(l2Tl). Curves show the runs of a^ t (Q) of 
Zhu et al. (2010a,b; green, short-dashed), a^.(Q) of Kratter et al. 
(2008; red, dotted), equation (|27l) . and aj)J(Q) of Martin &; Lubow 
(2011; long-dashed), equation (1261) . The latter is displayed for 
two values of the background viscosity: a m = 10 —4 (black) and 
otm = 10 -2 (blue.) 

One can see that reduces to a£ t if one uses 0.14 and 
1.3 instead of a m and Q C rit- Also, a^ t saturates at a 
constant value ss 0.1 for Q < 1, unlike which can be 
arbitrarily large for small Q. 

These a g t (Q) prescriptions are illustrated in Figure [U 
We summarize and analyze their properties in (© 

4. CONSTANT rn DISKS. 

One is often interested in knowing the structure of a 
steady-state accretion disk, which necessarily has M in¬ 
dependent of r. This limiting case provides a nice point 
of comparison between the different types of viscosity 
prescriptions and will be used in (© 

If M is independent of r (and, correspondingly, l = 
fir 2 ) and there are no external torques acting on the disk, 
then equation © naturally yields Fj = Ml (Rafikov 
2013). We now illustrate the difference in approaches 
between the various closure philosophies when calculat¬ 
ing the constant M disk structure. 

4.1. Constant M disk: Q = Q o closure 

When we explicitly demand a constant M disk to be 
marginally unstable with the condition (1141) , the equation 
(EUl) implies a relation between E, T, r in the form 

ipF J:gt (r, E) + F Jtm (r, E, T) = Ml(r). (28) 

The second algebraic relation between these variables is 
provided by equations ©-H3. thus fully specifying the 
behavior of E(r) and T(r). Rafikov (2009) and Clarke 
(2009) used this approach to understand the properties 
of gravitoturbulent disks with constant M. 
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4.2. Constant M disk: explicit a g t(Q) closure 

When one uses an explicit a g t(Q) prescription, the ap¬ 
proach to determining disk structure is different from 
that in H44.ll Assuming a non-zero background viscos¬ 
ity, equation I© now yields for the constant M disk 

37 tc 2 E [agt(Q) + a m ] = Mfi(r). (29) 

Here both c s and Q explicitly depend on temperature. 
Thus, agt(Q) is a function of T in this approach, unlike 
the case considered in H3 3.ll 
Again, equations ©-(HU) provide a second relation be¬ 
tween E, T, r, allowing one to fully determine the con¬ 
stant M disk structure (see e.g. Zhu et al. 2010a,b; Mar¬ 
tin & Lubow 2011). 


5. COMPARISON OF DIFFERENT O gt PRESCRIPTIONS 

We now compare the performance of different prescrip¬ 
tions for a gt described in 13 3.1113 3^1 To that effect we 
construct steady state (constant M) models of gravito- 
turbulent protoplanetary disks using the results of (© 
We then compare the behaviors of various disk charac¬ 
teristics obtained with different a g t recipes. 

Our calculations adopt the opacity behavior CEU-CE3 
typical for low temperatures T < 200 K (Bell & Lin 
1994). For simplicity we will assume this behavior to ex¬ 
tend also to higher temperatures; this should not be a 
problem as out present goal is to explore the differences 
between the various a g t prescriptions. We take the mass 
of the central star to be M* = Mg and set the threshold 
for gravitational stability at Qq = 1.5. In this and subse¬ 
quent calculations we assume that fragmentation ensues 
when a exceeds a cr it = 0.1 ( rathe r than unity, as e.g. 
is assumed in equations (lAH)-(lX3l)h Numerical studies 
suggest that a lower value of a cr i t is more realistic (Rice 
et al. 2005; Paardekooper et al. 2011). 

Rafikov (2009) showed that gravitoturbulent state re¬ 
sults in a set of fiducial values of various disk variables: 
surface density E f , midplane temperature and sound 
speed Tf and c s j, mass accretion rate Mf, and angu¬ 
lar frequency f If. They are chosen such that at the ra¬ 
dius corresponding to the angular frequency S2/ a steady 
gravitoturbulent disk with the mass accretion rate equal 
to Mf simultaneously (1) has midplane temperature and 
sound speed equal to Tf and c s j, (2) fragments, and (3) 
has optical depth r = 1. In Appendix [Al we provide ex¬ 
plicit expressions and numerical estimates for these vari¬ 
ables for the case of interest to us (opacity n oc T 2 ). 
These fiducial quantities provide characteristic values of 
the important parameters of gravitoturbulent disk and 
will be used in our subsequent comparisons. 

For a given mass of a central object M* characteristic 
angular frequency f If determines a fiducial distance r/ 
according to the formula 


r f 



130 AU 


( M+ \ 1/3 

\m q J ’ 


(30) 


where we took the numerical value of flf from equation 
(IA2I) in Appendix [A] This is the radius beyond which 
an optically thick disk with opacity re oc T 2 inevitably 
fragments (Matzner & Levin 2005; Rafikov 2009; Clarke 
2009). 


i o/n, 



Fig. 2.— Comparison of gravitoturbulent disk models computed 
using Q = Q o closure (black solid curve) and several different 
agt(Q) prescriptions: of Zhu et al. (2010a,b; green, short- 

dashed), q^J(Q) of Martin & Lubow (2011; blue, long-dashed), and 
a^.(Q) of Kratter et al. (2008; red, dotted). Plotted are: (a) opti¬ 
cal depth, (b) Toomre Q, (c) a-parameter, ( d) surface density, (e) 
midplane temperature. Opacity in the form (1121) - m is assumed. 
Background viscosity is high, a m = I O' 2 . A high-M case is con¬ 
sidered ( M = 10 Mf «7x 10 -5 Mg yr -1 ) with the disk remaining 
optically thick all the way out to Kb 150 AU. Beyond this radius it 
inevitably fragments as agt exceeds critical value a cr it = 0.1. 


In all our models we assume that the disk is immersed 
in a uniform radiation field with the radially-independent 
temperature T; rr = Tf ca 12 K. We have checked that all 
our conclusions remain the same for other values of Tj rr , 
in particular for the self-luminous disks with Xj rr = 0. 

In Figure [5] we show the radial profiles of the midplane 
temperature T, surface density E, optical depth r, a- 
parameter and Toomre Q for a gravitoturbulent disk that 
has a high mass accretion rate, M = 10 Mf w 7 x 10~ 5 
Mq yr -1 (leaving aside the question of how realistic such 
M is). We also use a relatively high value of the back¬ 
ground viscosity, a m = 10 , which is often adopted in 

the literature (Zhu etal 2009a; Martin & Lubow 2011, 
2013). Different curves correspond to a g t prescriptions 
given by equations (HU), CSD-KSZD, as indicated in the 
Figure. 

Shaded region outside of ~ 130 AU corresponds to a 
fragmenting part of the disk where a > a cr it, and the 
behavior of disk parameters in this region is irrelevant 
(many curves are very discrepant there). As expected 
(Rafikov 2009), at the fragmentation edge our disk model 
has E>E/,T>T/, and is optically thick (r ~ 10). 

One can see that all a g t prescriptions reliably repro¬ 
duce the radius at which the transition from the grav¬ 
itationally stable to the gravitoturbulent state occurs 
(around the vertical Q = Q o line). Interior to this radius 
background viscosity dominates, a ~ a m , and different 
curves overlap. But outside this radius, in the region 
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Fig. 3.— Same as Figure |21 but now for the low background 
viscosity a m = 10 —4 . One can see significant differences between 
the results obtained using a^(Q) recipe and all other a g t prescrip¬ 
tions, which tend to generally agree with each other. 


Fig. 4.— Same as Figure l3l (same low a m = 10 —4 ) but now 
for the low-M case, M = 0.1 Mf ~ 7 X 10 -7 Mq yr _1 . Disk 
becomes optically thin around 10 2 AU and remains stable against 
fragmentation out to > 10 3 AU. 


where a > a m , some differences between a gt prescrip¬ 
tions start emerging. In particular, at the fragmenta¬ 
tion edge (boundary of the shaded region) ct^j. computed 
through Q = Q q closure is different from of Martin 

6 Lubow (2011) by about a factor of 2. Similar differ¬ 
ences at the same r ss 130 AU are seen in the values of 
E and r. Other explicit a g t prescriptions show better 
agreement (within tens of per cent) for E, a, and r with 
the results of Q = Qo closure. On the other hand, the 
behavior of the midplane temperature (Figure [2^) shows 
good agreement between different a g t prescriptions. As 
expected, T never drops below about 10 K as it is lim¬ 
ited by Ti rr far from the star; however, for this disk model 
this is true only beyond the fragmentation edge so that 
in practice irradiation with the adopted T; rr is irrelevant 
here. 

To investigate these differences further we recomputed 
the disk structure for the same high mass accretion rate 

7 x 10 -5 M 0 yr" 1 but decreasing the background vis¬ 
cosity by two orders of magnitude, to a m = 10 -4 . This 
lower value of a m may be more realistic for the cold 
protoplanetary disks, where the MRI is weakened by the 
non-ideal effects such as ambipolar diffusion (Bai & Stone 
2011 ), or may not be operating at all, e.g. in a dead zone 
(Gammie 1996; Fleming & Stone 2003). 

Figure [3] shows the results of this calculation. Its in¬ 
spection reveals good quantitative agreement between 
the models computed via Q = Qo closure (black solid), 
which rely on given by the equation (fTSl) . and those 
using explicit viscosity in the form proposed by Kratter 
et al. (2008), a^(Q), and Zhu et al. (2010a,b), a^ t (Q). 
They typically agree with several tens of per cent for all 
variables in the whole gravitoturbulent region of the disk, 
between the fragmentation edge and the Q = Qo line. 


It is also obvious that the use of a^(Q) suggested by 
Martin & Lubow (2011) and a m this low leads to rather 
poor quantitative agreement with all other prescriptions: 
E, a and r show a discrepancy of more than an order of 
magnitude in the outer disk, at the edge of the frag¬ 
mentation zone. The value of Q at this radius plunges 
to ~ 0.2 for a^(Q), in contrast to Q « Qo maintained 
by all other prescriptions in agreement with simulations. 
Lower value of a m used in this calculation pushes grav¬ 
itoturbulent zone down to ~ 30 AU, accentuating the 
deviations which were already visible at higher a m , with 
less extended gravitoturbulent region (see Figure 0. 

This general situation holds for other disk models. For 
example, in Figure [I] we compare different a gt recipes in 
a low-M disk with M = 0.1 Mf = 7 x 10 -7 Mq yr _1 , 
while keeping a m = 10 -4 . This disk is optically thin in 
its outer parts, with r = 1 transition happening interior 
to rf. In this case, as shown in Rafikov (2009), fragmen¬ 
tation can be pushed out to very large radii and the disk 
can extend beyond 10 3 AU in the optically thin, gravito¬ 
turbulent regime. Moreover, external irradiation is very 
important in this case: 12 K sets the temperature 

floor outside 100 AU, in the part of the disk which is 
stable against fragmentation, see Figure [4jj. 

Again, with such low a m disk models using a^(Q) 
suggested by Martin & Lubow (2011) are considerably 
different from all other models, which tend to agree with 
each other. For E the discrepancy can be as high as 
an order of magnitude. From this we can conclude that 
a gt(Q) prescription underperforms at low values a m . We 
discuss the reasons for this next. 

6. DISCUSSION. 
















































Our interpretation of the behavior found in models us¬ 
ing a^l(Q) is that it is due to the rather gradual depen¬ 
dence of the former on Q. Indeed, Figure [1] shows that 
when a g t is between a m and the critical value for frag¬ 
mentation (a cr it = 0.1 in our case) a^(Q) grows with 
decreasing Q much slower than either a£ t (Q) or t {Q). 
This is especially obvious in the a m = 10 ~ 4 case, when 
Q has to go down to Q « 0.1 for a to reach a cr i t . This 
explains a significant drop in Q towards the fragmen¬ 
tation edge in Figures [ 3}3 and |4jr. Thus, we conclude 
that a^.(Q) does not properly capture the main feature 
of the gravitoturbulent disk - almost constant and close 
to unity value of the Toomre Q. The explicit dependence 
of a^(Q) on a m may be another feature that causes it 
to underpeform compared to a^ t (Q) and a^ t {Q). 

At the same time, we must emphasize that we are 
not aware of the existing gravitoturbulent disk calcu¬ 
lations using low a m = 10 -4 — most of them assume 
a m = 1CT 2 . Thus, Figures [3] and [4] are not providing 
comparison with any published results, but merely urg¬ 
ing caution in choosing a gt recipe when working with low 
a m . And at higher a m = 10 ~ 2 calculations using a^.(Q) 
can still be acceptable (see Figured]), depending on the 
required accuracy. 

Explicit prescriptions suggested in Kratter et al. (2008) 
and Zhu et al. (2010a,b) (a^.(Q) and a^ t (Q)) do rather 
well at reproducing the Q — Q 0 calculation using ajf t . 
Based on this we expect that effectively any explicit 
a gt (Q) scheme that exhibits very sharp rise in the vicin¬ 
ity of Qo as Q decreases, should be suitable for practical 
calculations of the gravitoturbulent disk properties. In 
fact, the sharper is the better and something as sim¬ 
ple as a gt (Q) = exp[(Q 0 — Q)/&Q} with ^ 10~ 2 will 
maintain Q « Q 0 in the gravitoturbulent state quite well 
(conditional upon a < a C rit)- 

Nevertheless, we strongly believe that there are signif¬ 
icant benefits to using the Q = Q o approach outlined 
in 33 3.11 First, the only assumption that it employs 
is that the disk is able to maintain itself in a state of 
marginal gravitational stability Q = Qo- This stipula¬ 
tion is in agreement with the results of direct numerical 
simulations (Cossins et al. 2009, 2010). Explicit a gt (Q) 
prescriptions also make this assumption, even though im¬ 
plicitly. But in addition they have to postulate some 
actual functional form of the a g t (Q) dependence, which 
does not have physical justification and is always intro¬ 
duced in an ad hoc fashion. This makes the explicit 
a gt (Q) approach rather arbitrary, which may result in 
certain problems, as we have demonstrated in f[5] 

Second, deep in the gravitoturbulent regime, when one 
can neglect the background viscosity a m compared to the 
gravitoturbulent contribution the explicit Q = Qo 
approach results in a self-contained equation m for the 
surface density evolution. This is allowed by the fact that 
the expression m already implicitly accounts for the 
thermal balance of the disk , even in presence of external 
irradiation. Disk temperature is uniquely defined in this 
case by the relation (1151) after the disk surface density 
has been solved for. This is not the case for the explicit 
a gt (Q) recipes, which always need to explicitly relate T 
to E via the equation of thermal balance (even deep in 
the gravitoturbulent state, when a gt a m ), because 



r, AU 

Fig. 5.— Behavior of a at the transition between the gravi¬ 
tationally stable and unstable disk regions. Curves of a (solid), 
a^. (dashed), a m (dotted) are shown as functions of r (lower axis) 
and Q (upper axis) for a m = 10 —2 (black) and a m = 10“ 4 (red). 
In both cases M = 10 Mf ~ 7 x 10 -5 Mq yr _1 , and these results 
correspond to Figures [2] and [3] One can see that in the gravitation¬ 
ally stable part of the disk (at smaller r and larger Q) rapidly 
becomes subdominant compared to the background visc osity am 
and a —> a m , thus justifying our simple prescription m . 


both enter the definition CD of the Toomre Q. Thus, our 
Q = Qo prescription allows a more efficient numerical 
exploration of the gravitoturbulent disk structure. 

Third, when the background viscosity cannot be ne¬ 
glected compared to a g t, the explicit Q = Qo prescrip¬ 
tion (fl8l) does not make the task of computing the disk 
properties more complicated than the explicit a g t(Q) 
prescriptions, see H3 3.1 3.1.11 Dealing with this case us¬ 
ing our prescription (1241) , which results in a non-zero a g t 
even in the gravitationally stable part of the disk, is not 
a problem. 

Indeed, in Figure[5]we demonstrate that a g t computed 
in our disk models via equation CHD rapidly becomes sub¬ 
dominant as the disk transitions into the regime dom¬ 
inated by the background viscosity a m . Thus, even 
though in general a^j. is not precisely zero 4 in this part 
of the disk, it still does not affect the disk properties. 
The only place where our approach may be quantita¬ 
tively less accurate is around the transition ct gt ~ a m . 
However, this is the location where all a g t prescriptions 
have some degree of arbitrariness. 

Fourth, Q = Qo prescription allows a much more trans¬ 
parent analytical derivation of certain gravitoturbulent 
disk characteristics. We demonstrate this next when we 
compare the derivations of the snow line location ( 1 16 6 . 11 ) 
and of the dead zone edge ( 36 6 . 21 ) in gravitoturbulent 
disks using the two classes of a g t prescriptions. 

4 Graviturbulent contribution does vanish inside of 25 A U in 
Figure [5l for a m = 10" 2 because of our use of b in the form (f23l 
and the fact that Tq drops below Ti rr inside this radius. 
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6 .1. Snow line location 

Snow (or ice) lines for different volatile materials (H 2 O, 
CO, etc.) arise in protoplanetary disks at the locations 
where these materials undergo a transition from solid to 
gaseous phase. It is conventional to identify snow lines 
with the disk radii where the midplane temperature T is 
equal to some characteristic value T snow for sublimation 
of a particular material at the typical midplane pressure. 
Following common wisdom we will adopt this as a work¬ 
ing definition of a snow line. 

Martin & Livio (2013) derived an analytical expres¬ 
sion for the snow line radius R sn ow in a gravitoturbulent, 
optically thick, self-luminous protoplanetary disk with 
constant M. Here we compare their calculation with 
the derivation of R sn ow using our favored a g t recipe (1181) 
based on the condition Q = Qo- To facilitate compari¬ 
son we employ a common set of assumptions: ( 1 ) neglect 
background viscosity a m and external irradiation (i.e. 
Tin = 0 ), ( 2 ) describe optically thick radiation transfer 
by 5 

/O) = (31) 

in place of our equation m, with t = (1/2)kE (instead 
of our definition m and (3) take opacity at the icy 
grain sublimation radius to be given by equation m 
with /3 = —0.01 and kq = 3 (in appropriate CGS units). 


6.1.1. I? S now via the Q = Qo prescription 

In our approach we first express E via R snow and T snow 
by setting Tq = T snow in equation (1151) . Then we plug 
this E(i?. sn 0 w) into the expression (fl 8 l) for a gt . Finally, 
inserting all this into equation (1281) . written in the form 
37rc 2 Ea gt = Mfl with all variables expressed through 
Rsnow and T sn ow, and resolving it for R sn ow we get 


= 9 kqAI 3 /- / Gk B \ ■ 

snow 1287T 2 aQo V f* ) 

/ ■ \ 2 /9 

M \ 


-\ 2 /9 


3 5.78 AU Q 0 


-2/9 


M \ 


1/3 


M 


e) 


T sr 


lO” 8 M 0 yr _1 
-0.78 


170 K 


Note that this derivation does not use thermodynamic 
balance condition m as it has already been used in de¬ 
riving the expression (fT51) for ajf t . 


6.1.2. R snow via the explicit a gt (Q) 

Now we outline the steps used in Martin & Lubow 
(2013) to derive i? snow . They adopted an explicit 
a g t(Q) = «oexp(— Q 4 ) with 6 ao = 10” 2 . This expres¬ 
sion is similar to that in Zhu et al. (2009b; see our equa¬ 
tion (12511 1 in spirit, but not in detail, as the maximum 
value of a g t(Q) is modulated by a small factor ao in the 
Martin & Livio (2013) case. Note that this anzatz would 

5 This is what equations (4) &: (5) of Martin & Livio (2013) 
effectively reduce to. 

6 Martin & Livio (2013) used a m rather than ao to denote a g t(0) 

in their equation (9) for a g t(Q), but this notation would cause 

confusion with the background viscosity in our paper. 


never allow a gravitationally unstable disk to fragment 
unless the critical value of a g t for fragmentation a cr it 
were very low, below ao- 

With this prescription the constant M assumption, 
embodied in the equation (1291) allowed Martin & Livio 
(2013) to relate Q at the snow line to the disk midplane 
temperature (assumed equal to T snow ) and R snow . Be¬ 
cause of the adopted form of a gt (Q) this relation de¬ 
pended on ao and involved a non-elementary function 
(Lambert function), complicating the analysis. The lat¬ 
ter problem was overcome by using an asymptotic re¬ 
lation for the Lambert function, valid for relatively low 
values of M < 10” 7 Mq yr” 1 . 

Additional (algebraic) relation between T snow , Rsno w, 
and E at the snow line comes from the thermodynamical 
relation |9|) with Fj = MVtr 2 . Eliminating E(i? sn 0 w) 
with the aid of these two relations Martin & Lubow 
(2013) obtained a scaling for i? snow in terms of powers 
of T snow , M, and the central star mass M, given by the 
equation (19) of their paper. 

6.1.3. Comparison of approaches 

Comparing the result for 7? snow in Martin & Lubow 
(2013) with equation (l32l) shows that the two are essen¬ 
tially identical, with all the scalings being the same and 
the numerical estimates different at 1% level for Q 0 = 1. 
However, it is clear that our derivation of the expression 
(1321) is much more straightforward and flexible. 

First, it involves only elementary functions. Second, 
it does not involve ad hoc factors, such as ao used by 
Martin & Lubow (2013) in their definition of a g t(Q), 
which are confusing and do not appear in the final result 
anyway. Third, our result (15^1) is clearly valid regardless 
of the value of M , while the asymptotic representation 
of the Lambert function used by Martin & Lubow (2013) 
works only for relatively low M < 10”' Mq yr . 

Analogous problems would arise if one were to use 
other explicit a g t(Q) prescriptions for analytical calcula¬ 
tions. This provides strong motivation for adopting the 
Q = Qo approach advocated in this work rather than the 
explicit a g t(Q) recipes, when studying gravitoturbulent 
disks. 


6.2. Outer edge of the dead zone 


In a very similar vein we can estimate the location of 
the outer edge of a dead zone (Gammie 1996) in a grav¬ 
itoturbulent protoplanetary disk. This edge should be 
located at the radius I?d, where the disk surface density 
is twice the surface density of the active layer E a a: 100 
g cm” 2 , down to which external ionizing radiation can 
penetrate into the disk and keep it fully ionized. 

For this estimate the equation m of the Q = Qo 
approach allows one to express the midplane temperature 
through i?d and E a , while cv]f t is already a function of 
the two, see equation (fl 8 l) . Then, repeating the steps in 
H6 6.1 6.lTTl and eliminating T(i?d) one obtains 


R d = (GM) 1/3 


; 21 AU Q 0 


4-/3 


9 kqM (k B /p) 

128tt (T eI” 2/3 (7rGOo) 2(4_/3) 

\ 1/9 


1/(15—3/3) 


4/9 


M 


lO” 8 M 0 yr 


-1 
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(M _\ 1/3 ( £ a V 1/3 

V M 0/ V 102 S cm_2 / 


(33) 


In making the estimate we again adopted the opacity in 
the form m , which is perfectly reasonable since Rd > 
-Rsnow and ice grains dominate opacity. 

It is easy to see that deriving Rd via some explicit 
a gt (Q) prescription would again require going through 
the rather contrived procedure described in H6 6.1 6.1.21 
involving non-elementary functions, their asymptotic ex¬ 
pansions, and so on. Thus, similar to N6 6.1 6.1^1 we 
can again argue that the Q = Qo approach provides the 
shortest and clearest path to finding Rd. 


7. CONCLUSIONS. 

One of the main outcomes of this work is the develop¬ 
ment of a self-consistent analytical formalism for the evo¬ 
lution of a gravitoturbulent accretion disk, including the 
possibility of external irradiation. It explicitly enforces 
Toomre Q in the disk to stay close to the value Qo needed 
for marginal gravitational stability. Implicit numerical 
implementations of this approach do exist (e.g. Terquem 
2008; Clarke 2009; Zhu et al. 2009a) but a complete an¬ 
alytical formalism was lacking until now. Our work fills 
this gap with the prescriptions described in 1 13 3.11 in 
particular the explicit expression (11811 for the gravitotur¬ 
bulent viscosity a^ t , which depends only on r and E and 
fully accounts for the thermodynamic equilibrium of the 
disk. Our formalism includes the possibility of the non¬ 
zero background viscosity a m ( 113 3.1 3.1.1|) provided by 
sources other than the gravitoturbulence (e.g. MRI), and 
external irradiation of the disk. This work thus general¬ 
izes the models of steady, constant M gravitoturbulent 
disks explored in Rafikov (2009). 

In fully gravitoturbulent state, when a m we 

derive an explicit equation m for the surface density 
evolution, which is fully self-contained — it involves only 
r and E — and is non-linear in general. This equation 


provides a powerful tool for exploring gravitoturbulent 
disk evolution and is thus an important result of our 
work. 

We then contrasted our Q = Qo approach with the cal¬ 
culations specifying an explicit dependence of the gravi¬ 
toturbulent viscosity ct g t on the Toomre Q, which is often 
done in the literature. We clearly demonstrated that our 
Q = Qo approach provides a much faster and more nat¬ 
ural route to deriving analytical expressions for various 
important characteristics of the gravitoturbulent disks, 
such as the locations of the snow lines (j j6 6.11) and the 
edge of the dead zone (16 6.21) . We also compared disk 
models computed using different a g t recipes ($5]) and 
found that some explicit a g t(Q) prescriptions are able to 
reproduce the results obtained using Q = Qo approach 
reasonably well. This is the case only for those prescrip¬ 
tions that cause ct g t(Q) to increase very steeply as Q de¬ 
creases in the vicinity of Qo (e.g. Kratter et al. 2008; Zhu 
et al. 2010a,b). Certain disagreement may arise if this re¬ 
quirement is violated (Martin & Lubow 2010), especially 
when the background (non-gravitoturbulent) viscosity in 
the disk is low (see our discussion in IQ. Thus, although 
it is often stated in the literature that the precise form 
of the a g t(Q) dependence used for building gravitoturbu¬ 
lent disk models is unimportant, we find this to be only 
partly true. 

Finally, we emphasize that our formalism has the min¬ 
imal number of imposed assumptions — namely, that 
Q = Qo is approximately maintained in the gravitotur¬ 
bulent state, in agreement with simulations. It thus pro¬ 
vides a natural and robust pathway to building efficient 
time-dependent models of gravitoturbulent disks around 
objects such as young stars, quasars, and Pop III stars. 


I am grateful to Zhaohuan Zhu for useful comments. 
The financial support for this work is provided by the 
NSF grant AST-1409524. 
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APPENDIX 


CHARACTERISTIC DISK PARAMETERS 


For our adopted opacity behavior m-m, assuming that fragmentation happens at a cr it = 1 , the fiducial values 
of the surface density £/, midplane temperature and sound speed Tf and c s j, mass accretion rate Mf 1 and angular 
frequency f If are 


E/ = 


c s,f — 


8 a 


M \ -7 


9 ttGQq J \kB 
9 nGQo f [i 


1/15 


14 gem 2 , Tf = 


/ 9 7tGQ 0 n 

\kB 

\8 aul , 

/i 


1/15 


12 K, 


8 ctkq \ ks 

8 

Mf = 3n 


1/15 


-^(ttGQq ) 4 (A 


0.22 km s 1 , 

-1/5 


T = 


8 a 

9 Kq 


nGQo 




1/3 


1.4 x 10" iU s“\ 


7.2 x 10 -6 M 0 yr -1 . 


(Al) 

(A2) 

(A3) 


The numerical estimates are for Qq = 1. More general expressions for these variables for different opacity behavior 
(/3 / 2) and a cr i t , as well as the coefficient in the expression (fl8l) for a gt , can be found in Rafikov (2009). 















